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A Summary of the Space-Time Conservation Element and Solution 

Element (CESE) Method 


Xiao-Yen Wang 

National Aeronautics and Space Administration 
Glenn Research Center 
Cleveland, Ohio 44135 


Summary 

The space-time Conservation Element and Solution Element (CESE) method for solving conservation 
laws is examined for its development motivation and design requirements. The characteristics of the 
resulting scheme are discussed. The discretization of the Euler equations is presented to show readers how 
to construct a scheme based on the CESE method. The differences and similarities between the CESE 
method and other traditional methods are discussed. The strengths and weaknesses of the method are also 
addressed. 

1.0 Introduction 

The space-time Conservation Element and Solution Element (CESE) method was introduced by 
S.-C. Chang (Ref. 1) in 1995 as shown in the Figure 1 timeline. In 1999, the two-dimensional (2D) 
triangular and three-dimensional (3D) tetrahedral unstructured Euler solvers were published (Refs. 2 and 
3), which was followed by the extension to the 2D quadrilateral and 3D hexahedral meshes by Zhang 
(Ref. 4) in 2002. The Courant number (i.e., Courant-Friedrichs-Lewy (CEL) number) insensitive scheme 
was presented in 2002 (Ref. 5) and 2003 (Ref. 6) and gave users a better way to minimize the numerical 
dissipation. In the same year, Huynh published a paper to address the similarities and differences between 
CESE and other methods (Ref. 7) and gave an insightful view of the CESE method. Local time stepping 
was introduced in 2005 (Ref. 8) for time-accurate Euler solutions, whereby different timescales are used 
in the computational domain and the time step at each mesh point can vary during the time-marching to 
achieve higher CFL numbers at all locations and all the time. This results in better accuracy (less 
dissipation) and less computational time. The 2D/3D Navier-Stokes solver was published by C.-L. Chang 
in 2006 (Ref. 9) and 2007 (Ref. 10). In 2010, the high-order (fourth-order) CESE scheme was introduced 
by S.-C. Chang (Ref. 1 1) and was extended to the 2D Euler equations in 2013 (Ref. 12) and 2014 
(Ref. 13). The strategy for handling meshes with a high aspect ratio (a viscous mesh near a solid wall) is 
described in Reference 14. 

The CESE method can be categorized as a cell-centered finite volume method. It aims to solve 
unsteady flow problems. The core ideas that motivated the development of the CESE method are 
(1) enforcing flux conservation over a discretized space-time domain, (2) constructing a nondissipative 
scheme as a baseline to fully control the numerical dissipation, (3) avoiding the Riemann solver and 
dimensional splitting by using a staggered mesh in space and time while maintaining robust handling of 
shock waves and discontinuities, (4) using a compact stencil (use immediate neighboring cells and/or 
points), and (5) using unstructured meshes for ease of dealing complex geometries. 
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S.-C. Chang 
(Ref. 1 ): layout for 
non-dissipative 
schemes and 
space-time 
formulations 


Zhang, Yu, and 
S.-C. Chang 
(Ref. 4): quadrilateral 
and hexahedral mesh 
CESE schemes 
S.-C. Chang (Refs. 5 
and 6): Courant number 
insensitive schemes 


Huynh (Ref. 7): 
address trade off 
between CESE 
and upwind 
schemes 


C.-L. Chang (Refs. 9 
and 10): edge-base 
derivatives approaches, 
triangular/tetrahedral 
viscous simulations 



Wang and S.-C. Chang 
(Refs. 2 and 3): 
staggered triangular/ 
tetrahedral mesh CESE 
schemes 


S.-C. Chang, Wu, 
Yang, and Wang 
(Ref. 8): time-accurate 
local stepping for 
one-dimensional 
flows 


S.-C. Chang 
(Ref. 11): high-order 
CESE methods for 
general conservation 
laws 


C.-L. Chang, 
Venkatachari, 
and Cheng 
(Ref. 12): local 
time stepping 
and high-order 
CESE method 


Figure 1. — Timeline of CESE method development. 


The CESE method has been used to solve computational aeroacoustic (CAA) benchmark problems 
involving fan noise and jet noise (Refs. 15 to 20). All problems were successfully solved as blind tests; 
numerical solutions were compared with analytical solutions and showed good agreement, demonstrating 
the ability of the CESE method to solve complicated flow physics (e.g., shocks + waves + eddies). In 
References 21 to 25, the capability of the CESE method for solving engineering problems including high- 
speed flows with shock waves are demonstrated. 

However, some issues remain, such as 

(1) The CESE solver becomes very dissipative when the CFL number is small («1); here, either a 
CFL-insensitive scheme might help (Refs. 5 and 6) or local time stepping for a time-accurate solution 
would be possible way to remedy this (Ref. 8). 

(2) A highly stretched viscous mesh near the wall might degrade the accuracy of the gradients of the 
marching variables. There are a few ways to get better accuracy of the gradients; these can be found in 
References 7 and 14. Recently, a higher order CESE scheme (fourth-order) has been developed (Refs. 1 1 
to 13), which might enhance accuracy for the viscous flow problems that are dominated by boundaries 
that are solid walls. 


2.0 Discretization 

In the following, the numerical scheme is constructed for the 2D and 3D Euler equations using the 
CESE method. The governing equations are introduced, then the control volume is described in 
comparison with the traditional methods, and the CESE 2D and 3D Euler solvers are presented. 
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2.1 The Governing Equations 

The 3 D Euler equations in conservation form are 


8u m + + cfm q 

dt dx dy dz 


m = 1,2, 3 , 4, 5 


( 1 ) 


where u m are the conservative variables and and g m are the x-, y-, and z-component inviscid fluxes 
that are defined as follows: 


u i= p, U2= pu, u 3= pv, 114= p w, u$= pe t (2) 

2 

e\ = pu, 62 = pu' + p, e 3 = p uv, 64 = puw, 65 = (pe, + p)u ( 3 ) 

f\ = pv, fi = pvu, f 3 = pv 2 +p,f. 4 = pvn; f 5 = (pe t +p)v ( 4 ) 

gi = pm, g 2 = pwu, g 3 = pwv, g 4 = pm 2 +p,g5= (pe, + p)w ( 5 ) 


Here p is density, (u,v,w) are velocity components,/) is pressure, and e, is the total energy per unit mass. 
Its integral form for a finite control volume V in space (x,y,z) is 


IL 


du m + dem_ + ^fm_ + dg„ 


>vy dt dx dy dz 

Defining h m = ( e m ,f m ,g m ) in (x,y,z), Equation (6) then becomes 

du 


dE = 0 


( 6 ) 


II, ,^r ir +\ib -~ h -) ir = 0 


( 7 ) 


with V = 


r d d 


^dx dy dz 

into surface integration: 


. Using Gauss’s theorem, the volume integration in Equation ( 7 ) can be converted 


[[ ^dV+h h m -ds=0 ( 8 ) 

JJf dt ■> s(P) 

where s(V) is the boundary of the finite volume V. 

Equation (8) can be discretized using different traditional numerical methods. The volume integration 
of the time derivative allows people to use a different discretization approach for the time derivatives, 
which may be desirable for saving computational time, obtaining higher accuracy for steady-state 
solutions, or maintaining flux conservation in both space and time for transient problems. 

Contrary to the standard approach outlined in Equations (6), ( 7 ), and (8), with the CESE method we 
define h m = (e m , f m ,g m ,u m ) in space and time (x,y,z,t), and Equation (6) becomes 

JJ [ (v./? m )dF = 0 (9) 
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f 


and Fis a finite control volume in space and time. Next, using Gauss’s 


where V 


_e_ _5_ _a_ 

Jdx ’ 8y ’ dz’ dt) 


theorem and converting the volume integration in Equation (9) into surface integration: 



h m ■ d? = 0 


( 10 ) 


The CESE method is about the discretization of Equation (10). Note that the integral form of the 
Euler equation in Equation (10) is the original form of conservation laws. Equations (1) and (7) can be 
derived from Equation (10) only if the flow is continuous. For the time derivative term, CESE does 
surface integration instead of volume integration, which results in no ambiguity for the discretization of 
the time derivative term. Also, the scheme is second-order accurate in both space and time. The flux 
conservation in both space and time is guaranteed. 

In the following section, the discretization of Equation (10) will be described using a staggered 
(CESE method) versus nonstaggered (traditional method) control volume. A 2D unstructured triangular 
mesh is used as an example to explain how to construct a scheme for Equation (10) based on the CESE 
method in comparison with traditional methods. 


2.2 Control Volume (Traditional Method Versus CESE Method) 

Figure 2 shows an unstructured triangular mesh as an example. Consider a triangle, ABDF, which has 
three neighboring triangles. Let points A, C, and E be the centroids of those three neighboring triangles. 
With the traditional method, each triangle will be the control volume where the flux is assumed to be 
conserved. The solution point is at the center of each triangle. Then an interface exists between 
neighboring control volumes. Each interface is associated with multiple solution points. Therefore, the 
uniqueness and smoothness of the interface flux has to be addressed. Most of the time, the interface flux 
is computed using the solution point on the upwind side, a mixture of both sides, or an average of both 
sides, which results in the various classes of upwind and central differencing schemes. 

In the CESE method, for the triangle ABDF, the control volume is the hexagon ABCDEF formed by 
connecting the three neighboring centroids to the vertices of the triangle ABDF. One of the neighboring 
triangles of ABDF is ADIF, and the corresponding control volume is the hexagon FGDHIJ. Note that 
there is no interface between the control volumes of two neighboring triangles, but the two control 
volumes overlapped on GDEF. The solution point will be the geometric center of the hexagon, not the 
geometric center of the triangle. Each triangle is associated with a hexagon and a solution point, which 
again is the centroid of the hexagon. 

The conservation element (CE) is the same as the traditional control volume but in both space and 
time. The solution element (SE) is the boundary of the CE for each solution point, also in space and time. 
For an example, the CE for point G is the hexagonal cylinder with the base of the hexagon ABCDEF, and 
the SE will be the boundary surfaces of the hexagonal cylinder. 




(a) (b) 

Figure 2. — Control volumes for discretization of Euler equations, (a) Traditional method, (b) Conservation Element 
and Solution Element (CESE) method. 
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Within the SE associated with a solution point (j,n ), the numerical approximation of 
h m = (e m ,f m ,g m ,u m ) is denoted as h m =(e m ,f m ,g m ,u m ) with 

u m (x, y,z,t; j, n) = (u m ) j + (u mx ) / (x-x / ) + (u my ) j (y -yj) + j (z - zj ) + (w mf )j(t-t n ) (11) 

e* (x, v, z,t;j, n) = (e m ) 7 + (e,„ x ) / (x-x / ) + (e„ lv )/ ( v - Vy ) + (e mz ) y (z - z 7 - ) + (e ml ) j(t-t n ) (12) 

(-^5 v 5 j y n) = (fm )_/ + ( fnx ) j if my) j ( V — Vj ) + (fmz) j (z ~ Zj) + (fmt) j(t ~ t ) (13) 

g* (x, V, z, t; j, n) = (g m ) j + (g„ u .) (x — Xy ) + (g mv ) j ( v - Vy ) + (g OTZ ) j (z. - Zj ) + (g mt )j(t~t ") (14) 


For the CE at (j,n), Equation (10) becomes 

f> h*„ ■ d5 = 0 (15) 

Js(CE) 

Details on solving Equation (15) will be given in Section 2.5, Three-Dimensional Formulation. Flere 
the difference between the traditional and CESE methods is summarized in Table I, using ABDF as an 
example. 


TABLE I.— COMPARISON BETWEEN TRADITIONAL (UPWIND) AND CONSERVATION ELEMENT 
AND SOLUTION ELEMENT (CESE) METHODS 
[Points A through G are illustrated in Figure 2.] 


Definition 

Upwind scheme 

CESE scheme 

Comments on CESE 

Control volume (CV) 

ABDF (triangle, 
nonoverlap) 

ABCDEF (hexagon, 
partially overlap) 

Inflated CV, but the same stencil 

Boundary of CV 

BD, DF. and FB (3) 

AB, BC, CD, DE, EF, and 
FA (6) 

Doubled number of boundaries 

Number of interfaces between 
CVs of two neighboring 
triangles 

Three 

Zero (no interface) 

No need to consider interface 
fluxes 

Number of solution points each 
boundary of CV is associated 
with 

Two 

One 

Solution point is on the boundary 

How flux at center of boundary 
is computed 

Taylor series expansion 
from left and right and 
then upwind flux 

Taylor series expansion 
from the solution point 

Uniquely defined 

Example 

For BD, use solution at 
either C or G or both 
based on upwind 

For BC (and CD), use the 
solution at C 

Sum fluxes through all boundaries 
to update solution at G 
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2.3 Time-Marching for the Primary Unknowns (w m )” +1/ With CESE Method 

Figure 3 shows the control volume associated with the triangle ABDF in space and time (x,y,t). Given 
the solution at points A, C, and E at time level n, the solution at point G’ at t = n+ 1/2 is updated. Each 
solution point at t = n is associated with two surfaces of the control volume along the time direction and 
one surface at a constant time level. The flux is computed at the geometric center of each surface from the 
solution point based on a Taylor series expansion. For the surface on the top (time level n+ 1/2), the flux is 
simply the solution at the geometric center of the hexagon multiplied by the surface area. Some details are 
given as follows. 

(1) Considering the neighboring point A(ji,n) 


Compute h m = ( e m ,f m ,u m ) at the centroid of surface A’B’BA (A’F’FA ) using the Taylor series 
expansion from the solution point A as shown in Equations (1 1) to (14), then multiply the surface 
area, As^i ( As / 1j2 ), and its normal vector in the form of n = (n x ,n v , 0) . For the surface of ABGF, 

where time is constant, again compute h, n =(em,Jm,Um) at the centroid of the surface ABGF that has 
the normal vector in the form of n - (0,0,1) and an area of A.syi o , and then we have the flux 

2 

R'i i=2>A T f nifty) j\,k ^ As'yiy + (fr/d/I.O ^ AVy]_() (16) 

k = 1 


(2) Repeat the same procedure for the second neighboring solution point C(/ 2 ,«) to get the flux (R'ji) 
through the surface of B’C’CB, C’D’DC, and CDGB. 

(3) Repeat the same procedure for the third neighboring point E(/ 3 ,«) to get the flux (Z^/V) through 
E’D’DE, E’F’FE, and EFGD. 


NASA/TM— 2015-218743 


6 


(4) At t = 72+1/2, the surface A’B’C’D’E’F’ whose area is denoted as A Sj, the flux is simply 
As 'j( u m )"j +ll ~ because the solution point is located at the geometric center of the hexagon 
ABCDEFG. 

(5) Summing all fluxes through the CE for solution point G’(/, 77+1/2) based on Equation (15) will 
give 


A Sj{u m )f 12 = -(/?;, +R] 2 +R%) (17) 

which gives the explicit time-marching scheme to update ( 2 . The stability condition of the 
scheme is CFL < 1. 


2.4 Update the Gradients 

A weighted-average limiter based on the gradient amplitude is used to compute u mx and u mv at G using 
the three neighboring solution points (A, C, and E) (Fig. 4). 

Three sets of gradients, such as ( iiml,iiml ), where / = 1,2,3, can be computed using the three planes 
GAC, GCE, and GEA using a finite-difference approach. Define q mJ = ^(u^ly + (Umly , and 
W»,1 = (q m ,2 < im,3) a ,w m ,2 = ( q m ,\qm,3) a ,w m ,3 = (q m ,\q m,2 )“, where a = 0, 1,2. The weighted-average limiter is 


(1) , (2) , (3) 

W m ,lUmx + W m .2U mx + . 3 ^ mx 

W m , 1 + W m ,2 + W;h,3 


- , U 


my 


(1) , (2) , (3) 

W w j 22 my + ^m,2^my + W m ,2Mmy 

W m , 1 + V + Wm,3 


( 18 ) 


Equation (18) is the original way of computing the gradients proposed by S.C. Chang (Ref. 2). In 
Reference 6, S.C. Chang describes a modified version of Equation (18) to provide better accuracy of 
gradients of the time-marching variables. 



Figure 4. — Stencil used to compute derivatives. 
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Figure 5. — Subconservation element related to three-dimensional tetrahedral mesh. 


2.5 Three-Dimensional Formulation 

Consider a tetrahedron ABCD in Figure 5, in which one of its four neighboring tetrahedrons ABCP is 
also plotted. The tetrahedron shares the face ABC. Points G and FI are the centroids of ABCD and ABCP, 
respectively. The polyhedron GABCFI is the subconservation element (sub-CE) associated with solution 
point G. 

Similarly, three additional sub-CEs associated with point G can be constructed by considering in turn 
the three tetrahedrons that share one of its other three surfaces with ABCD. The summation of the four 
sub-CEs will be the CE, the control volume where flux is assumed to be conserved. It includes the 
tetrahedron itself (ABCD) and one-fourth of each of its neighboring tetrahedrons. The solution point is 
located at the geometric center of the CE, which is not necessarily at the center of the tetrahedron. 

The control volume in space and time (x,y,z,t) cannot be plotted. The solution at point G’ (/',«+ 1/2) is 
updated using four neighboring solution points at time level n. Each solution point at t = n is associated 
with three volumes of the control volume along the time direction and one volume at constant time. The 
flux is computed at the geometric center of each volume from the solution point based on a Taylor series 
expansion. At the volume on the top (time level 77+1/2), the flux is simply the solution at the geometric 
center of the CE multiplied by its volume. The fluxes are then summed over the control volume based on 
Equation (15). Some details are given as follows. 

1 . Considering the neighboring point Yi(j\,n) 

Compute h m =(e m ,fm,gm,u m ) at the centroid of volume HABH’A’B’ (HACH’A’C’ and 

HBCH’B’C’) using a Taylor series expansion from the solution point H, then multiply the volume 

Avji ji , where k = 1,2,3 and the normal vector of the surface HAB (HAC, HCB) that is in the form of 

h = (n x ,n v ,n : ,0) . (H’,A’,B’,C’ is the extrusion of the point H,A,B,C in the time direction, not 

plotted in the figure.) For the volume GABCH, where time is constant (t = n), again compute 
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hn at the centroid of the volume GABCH that has the normal vector in the form of 

77 = ( 0 , 0 , 0 , 1 ) and volume of Av 7 -| i0 . W e have the flux 

3 

R‘n =2>k T fnJly “I" gmft r) / l . k * AVyl t k “1“ (Mm) / 1 .0 ^ Av / 1.0 (19) 

k = 1 


2. Repeat the same procedure for the rest of the volumes associated with the other three neighboring 
solution points to get the fluxes Rj 2 ,Rji, and Rj 4 . 

3. At t = n + 1/2, the volume of the CE is A v j , and the flux is simply (w m )" +1 2 • A vj because the 
solution point is located at the geometric center of the CE. 

By summing all fluxes through the CE for solution point G ’(/, n+1/2), Equation (15) yields 

A Vj(u m )f 12 = -(Rji + Rji +R’h + Rja) (20) 

which results in the explicit time marching scheme for (w w )y +1/2 . The gradients can be easily extended to 

3D formulation based on Equation (18) and will not be repeated here. 

The location where the flux is computed at the boundary of the CE is summarized in Table II to give 
readers a better idea how it is defined. 


TABLE II.— LOCATION OF THE FLUX IS COMPUTED IN ONE. TWO. AND THREE DIMENSIONS 


Dimension 

Boundary of conservation element (CE) along the time direction 

Boundary of CE at t = n 

3D 

Three surfaces, geometry center of each surface, time is at t = f + d7/4 

One volume, center of the volume, t=f 

2D 

Two lines, geometry center of each line, time is at 7=7" + d//4 

One surface, center of the surface, t = t" 

ID 

One point, at that point, time is at t = t" + dt/4 

One line, center of the line, t = t" 


2.6 CESE Nondissipative Scheme 


The schemes in Equation (17) for 2D and Equation (20) for 3D are used in most computations of 
Euler flow. With the weighted-average limiter, the numerical dissipation is adjusted from the 
nondissipative scheme with added numerical dissipation as needed. The nondissipative scheme is the 
foundation of the CESE scheme and guides the CESE scheme construction under a set of rigorous rules 
(Ref. 1). When solving real problems, how much dissipation is needed depends on the physics of the 
problem. In general, numerical dissipation has to be less than the physical viscosity in order to obtain the 
correct physical solution. Sometimes the physical viscosity is extremely small, and then it requires the 
scheme to be highly accurate yet stable. 

For the 2D case described in Section 2.2, Control Volume (Traditional Method Versus CESE 
Method), the conservation element in space is defined as the hexagon ABCDEF. It can be seen that the 
hexagon is formed by three quadrilaterals, ABGF, BCDG, and DGEF, that can be defined as three sub- 
CEs. Assuming Equation (15) is true for the three sub-CEs, there will be three equations for three 


unknowns (u m )’j +1/2 ,(u. 


mxj j 


\7?+l/2 1/ xW + 1/2 

) j , cind (u m y) j 


for each m. For the 3D case, each conservation element has 


four sub-CEs that yield four equations for four unknowns (u m ) 


n+l/ 2 

j ' 


/ \72+l/2 / 

y^mx ) j ? W 


my J j 


xH+1/2 i / \72+l/2 

) j ? and (u mz ) j 


for each m. The marching variable (u m )j' and its spatial derivatives are solved using the conservation 
laws for the nondissipative CESE scheme. 
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2.7 CESE Fourth-Order Scheme Compared to Traditional Methods 

The fourth-order CESE Euler solver has been developed recently. Validation and assessment is under 
way (Refs. 12 and 13). The stencil and control volume used in the fourth-order scheme is the same as that 
used in the second-order scheme. The time-marching scheme is the same as that of the second-order 
scheme, but it is used twice. First, the differentialized governing equation for the higher order derivatives 
(Uxx and Uxxx) is solved; then the original governing equation for u and u x is solved. Fligher order Taylor 
series expansion is used to compute the flux on the boundary of control volume. The stability condition is 
CFL < 1, the same as for the second-order scheme. The shock-capturing limiter is the same as that for the 
second-order scheme. Upwinding is not needed since the flux is uniquely computed at the interface. The 
traditional high-order finite-volume methods use a high-order polynomial to represent the solution inside 
each cell. They (1) divide each cell into the subgrid based on the order of accuracy, (2) decide the 
polynomial coefficients using the subgrid solution, and (3) reconstruct the flux on the interface of the 
control volume. Note this has to be carefully handled since the accuracy of the scheme is highly 
dependent on this. The higher the order of accuracy, the more restricted the stability conditions become; 
that is, the smaller the time step size will be. 

3.0 Discussion 

The CESE scheme has some concepts that are similar to those that were previously described in the 
literature. 

Space-time formulation was mentioned in Roe’s paper in 1983 (Ref. 26), but the importance of the 
space-time formulation for the scheme’s accuracy was not discussed. CESE emphasizes space-time flux 
conservation and applies it to the discretization scheme. Given that flux conservation is applied in space 
and time, there is no ambiguity on how to discretize the time derivative terms with the CESE method. 

A staggered mesh was used by Nessyahu and Tadmor (N-T) on a 2D quadrilateral mesh in 1990 (Ref. 
27); but it was not extended to a triangular mesh for 2D or any 3D formulation. The CESE method is 
consistent in its formulation for one, two (triangle), and three (tetrahedron) dimensions. Also, the N-T 
scheme does not have a nondissipative scheme. The CESE nondissipative scheme is the core of the 
scheme development. 

The weighted-average limiter was proposed by van Albada in 1982 (Ref. 28), but it is for a ID 
formula; the triangular mesh version of the corresponding limit was not developed. The CESE method 
has a coherent formulation of the weighted-average limiter from one, two, and three dimensions, which 
works robustly for capturing discontinuities, such as shock waves and contact surfaces. 

Staggered versus nonstaggered control volumes is the major difference between the CESE and 
traditional methods. Although a staggered control volume avoids the interface so there is no need of 
upwinding for the flux at the interface, its effect on the overall performance of the scheme compared with 
that of a nonstaggered control volume needs to be further evaluated using real engineering problems; that 
is, 3D unsteady viscous turbulent flows. The performance of the scheme includes the cost (computation 
time), accuracy, and robustness. These efforts are underway. 


4.0 Concluding Remarks 

The Conservation Element and Solution Element (CESE) method for solving conservation laws has 
been described and compared with other traditional methods. The numerical schemes for both two- and 
three-dimensional Euler equations were illustrated based on the CESE method. The major differences and 
similarities between the CESE method and other traditional methods have been highlighted and addressed 
from different perspectives. 
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